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Context. In stellar interiors, rotation is able to drive turbulent motions, and the related transport processes have a significant influence 
on the evolution of stars. Turbulent mixing in the radiative zones is currently taken into account in stellar evolution models through a 
set of diff'usion coefficients that are generally poorly constrained. 

Aims. We want to constrain the form of one of them, the radial diff'usion coefficient of chemical elements due to the turbulence driven 
by radial differential rotation, derived by Zahn (1974, 1992) on phenomenological grounds and largely used since. 
Methods. We performed local, direct numerical simulations of stably stratified homogeneous sheared turbulence using the Boussinesq 
approximation. The domain of low Peclet numbers found in stellar interiors is currently inaccessible to numerical simulations without 
approximation. It is explored here thanks to a suitable asymptotic form of the Boussinesq equations. The turbulent transport of a 
passive scalar is determined in statistical steady states. 

Results. We provide a first quantitative determination of the turbulent diff'usion coefficient and find that the form proposed by Zahn 
is in good agreement with the results of the numerical simulations. 
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4^ 1. Introduction 

I As many astrophysical measures are calibrated on stellar evolu- 
O tion theory, realistic stellar evolution models are crucial for as- 
^ trophysics. One of the major issues is the influence of the macro- 
^ scopic motions induced by the rotation of a star on its internal 
■ ■ structure (see the recent review of Maeder & Mey net 2012 ) . In 
particular, chemical mixing can significantly extend the lifetime 
of the stars by continuously providing fresh combustible matter 
to the core, thus feeding nuclear reactions for a longer time. 

A key point of these stellar evolution models is the radial tur- 
bulent mixing of chemical species induced by radial diff'erential 
rotation. According to Zahn's model ( |Zahn|1974| ), this turbulent 
mixing is described by a diff'usion coeflScient that reads as 
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where k is the thermal diff'usivity, Ricr the critical value of 
the Richardson number, Nj the Brunt- Vaisala frequency, and 
rsin^dD/dr the shear associated with radial diff'erential rota- 
tion. This expression is taken from Eq. (2.14) of |Zahn| ( [l992| ). 
While other versions of this difl'usion coefl&cient exist in the lit- 
erature, notably those accounting for additional effects like the 
//-gradient (Maeder & Meynet 1996;,Talon & Zahn|1997l|Mathis| 



et al.|20Q4| , Zahn's prescription is the basis of most implemen 



tations of the turbulent mixing driven by differential rotation in 
stellar evolution codes (|Endal & Sofia] 1978t [Pinsonneault et al 



1989^ "Heger et al.'2000' see a detailed discussion in Meyn et & 
Maeder 2000). Moreover, this mixing plays a significant role as 



compared to other transport processes, such as the meridional 
circulation combined with horizontal turbulence, especially for 



massive stars ( [Meynet & Maeder] 1997 ','2000 ; Denissenkov et al. 



1999| ) and low-mass stars with strong //-gradients ( jPalacios et al. 



2003 [ |2006| ). It thus has a direct impact on the confrontations 
of stellar evolution models with observations. While chemical 
abundances at the star surface have been widely used for such 
confrontati ons ( [Michaud & Zahn| 19981 [Brott et al.|2011t [Potter | 
etal. 12012] ), asteroseismology is now providing direct constraints 
on these processes, especially for subgiant and giant stars ([De-[ 
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[heuvels et al.|2012[[Eggenberger et al.|2012[ ). 

The derivation of Zahn's prescription is based on the hypoth- 
esis that turbulent flows are statistically stationary and that this is 
achieved when the mean shear flow is close to the linear marginal 
stability criterion. Following [Townsend] ([1958), Zahn modified 
Richardson's instability criterion to take the destabilizing eff'ect 
of thermal diff'usivity into account in such a way that marginal 
stability is controlled by RiPe ~ 1, where the Peclet number 
Pe compares the diff'usive time scale with the dynamical one (as 
later confirmed through numerical linear stability analyses, e.g. 
Lignieres et al.[p^99] ). Then, he assumed that, in a turbulent 
state, the relevant Peclet number in the expression RiPe ~ 1 is 
the turbulent one Per = — , defined with the turbulent scales of 
velocity and length (u and respectively). The diff'usion coef- 
ficient is then estimated through the relation D = u£/3, which 
finally leads to Eq. ([T]). 

The purpose of our present work is to use local 3D direct 
numerical simulations (DNS) to test Zahn's modeling of turbu- 
lent mixing. Previous numerical simulations performed in a geo- 
physical context have shown that a steady stably stratified homo- 
geneous sheared turbulence can be obtained by imposing both a 
uniform mean shear and a uniform stratification and by tuning 
the Richardson number to reach statistical stationarity ( J acobitz[ 
[et al.[[1997[). Such a flow is well adapted to studying chemical 



turbulent transport and enables one to relate the transport to the 
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as length, velocity, and temperature units (see Fig.[T]), their non- 
dimensional form reads as 



Fig. 1. Sketch of the flow configuration. 



local gradients of velocity and entropy as in Zahn's criterion. Im- 
posing the shear and the stratification amounts to an assumption 
of scale separation between the variations in the mean and tur- 
bulent flows. The main difi&culty of performing such numerical 
simulations in a stellar context is the very low Prandtl number of 
the stellar fluid. Indeed, this introduces a huge gap between the 
dynamical time scale and the diflTusive one, which forces use of a 
large number of numerical time steps to accurately compute the 
efiTect of thermal dififusion on a dynamical time. This would re- 
quire a prohibitive amount of computation time for DNS, where 
the whole range of scales of turbulence down to the dissipative 
scales ar e resolved. An attempt to test Zahn's prescription was 
made by |Bruggen & Hillebrandt| ( [2001| ), but in this case the sim- 
ulations relied on the numerical diflTusivity of the ZEUS code that 
is strongly dependent on the grid size (Fromang & Papaloizou 



[2007 ). Fortunately, an asymptotic form of the Boussinesq equa- 
tions in the limit of the small Peclet numbers exists, where the 
numerical constraint associated wit h the dynamical efiTect of the 
high thermal diflTusivity disappears ( |Lignieres|1999 ). 

In this paper, we perform DNS of 3D steady, stably stratified, 
homogeneous sheared turbulence at decreasing Peclet numbers 
and use both the Boussinesq equations and the so-called small- 
Peclet-number approximation (SPNA) to explore the asymptotic 
Pe ^ I regime. The mathematical model of the considered tur- 
bulent flow is presented in Sect. [2] followed by the numerical 
simulations (Sect. [3]) a nd the methods used to study the turbu- 
lent transport (Sect. [4|). Then, in Sect. [5] constraints on the dif- 
fusion coeflftcient obtained from our numerical simulations are 
described and compared to Zahn's prescription. The results are 
discussed in Sect.|6l 
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where u = uex + vey -\- we^ is the velocity, p the pressure deviation 
from hydrostatic equilibrium, the temperature deviation from 
the mean profile, and fy and forcing terms (see Sect. 2.2 ). 

These equations show three non-dimensional parameters: 
(i) the Richardson number Ri = (Nj/S)'^, noting S = dU/dz = 
AU/L the shear rate, A^t^ = Q;gAT/L the Brunt- Vaisala fre- 
quency, and a the thermal expansion coeflftcient, (ii) the "macro- 
scopic" Reynolds number Re = LAU/v = SL^/v, where v is the 
molecular viscosity, and (iii) the "macroscopic" Peclet number 
Pe = LAU/k = SL^/k. 

To explore the regime of very low Peclet numbers, we use 
the SPNA ( Lignieres 1999 ). The basic principle of this approx- 
imation is to Taylor-expand all variables up to first order in Pe. 
Thus, Eqs. ([3]) and ^ become 
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(6) 



with ifr = 6/Pe. The background thermal stratification remains 
fixed because temperature deviations become infinitely small as 
Pe approaches zero. There are two non-dimensional parameters 
left: the Reynolds number and the "Richardson-Peclet" num- 
ber RiPe, which compares the eflfect of the stratification modi- 
fied by the thermal diffusion (with a time scale of tb^/t^, where 
tb = and = L^/k are the buoyancy and diffusive time 
scales, respectively) to that of the shear. This number is known 
to control the linear stability of shear flows at low Peclet number 
( Lignieres et al.|1999 , and references herein). 

An expansion similar to the SPNA has been proposed for an 
unstable thermal stratification (Thual 1992). However, it does 
not allow the feedback of the convective turbulent motions on 
the stratification, as observed in the solar convective zone. 

The numerical code uses a Fourier colocation method in the 
horizontal directions, compact finite differences in the vertical, 
and a projection method plus Runge-Kutta for time advancing. 



2. Mathematical model of the flow 

To obtain results that are as generic as possible, we consider the 
simplest flow likely to produce stably stratified, homogeneous, 
sheared turbulence. Its characteristics are a uniform vertical 
shear of horizontal velocity and a uniform vertical temperature 
gradient, as shown in FigJT] The governing equations of the flow 



are presented in Sect. 2.1 and the boundary conditions along 



with the forcing and the initial conditions in Sect. 2.2 



2.1. Governing equations 

We use the Boussinesq equations in which density fluctuations 
are neglected, except in the buoyancy term. For the usually low 
Mach-number flows inside stars, this approximation is justified 
as long as the motion vertical length scale is smaller than the 
pressure scale height. Using L, AU, and AT = (r^ax - ^rmn)/2 



2.2. Boundary conditions, forcing, and initiai conditions 

Periodic boundary conditions apply in the horizontal directions, 
whereas in the vertical direction the fluid is bounded by two hor- 
izontal surfaces that are impenetrable (w = 0), stress-free in the 
span wise direction (d^u = 0), uniformly sheared in the stream- 
wise direction (d^v = S), and thermally controlled (6 = 0). 
The mean shear profile U(z) is imposed thanks to a body force 
fy = [U(z) - vey]/At applied at every time step At, where v is the 
horizontal average of the streamwise velocity. The mean verti- 
cal temperature profile is imposed in the same way. With such 
a set-up, turbulence is not rigorously homogeneous in the verti- 
cal direction. Nevertheless, as already found in Sc humacher &| 
|Eckhardt| ( |2000| ) using similar forcing and boundary conditions, 
non-homegeneity is limited to thin layers (~ 10% of the domain 
vertical extent) close to the upper and lower boundaries. 

As initial conditions, we construct a statistically homo- 
geneous and isotropic Gaussian random velocity fluctuation 
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Fig. 2. Evolution of the turbulent kinetic energy K for three different 
Richardson numbers at turbulent Peclet number Pe^ = 52. After a tran- 
sient regime, the turbulent kinetic energy is stationary for Ri = 0.124, 
increases for Ri = 0.05 and decreases for Ri = 0.2. 



field ( [Qrszag|[l969| ) with an energy spectrum proportional to 
k^e~^ /^o ^ where k is the norm of the wave vector and cor- 
responds to the maximum of the spectrum ( [Jacobitz et al.|1997] ). 
The initial temperature deviation 6 is set to zero. 

3. Numerical simulations 

For each value of the Peclet number, the Richardson number can 
be tuned to obtain statistical steady- state flows. In the regime 
where thermal diffusivity has no dynamical eff'ect, various stud- 
ies ( |Holt et al.|p^92l [Jacobitz et aLjpWT] ) had already shown 
there is a critical Richardson number such that the turbulent ki- 
netic energy increases if Ri < Ricr and decreases if Ri > Ri^v, as 
illustrated in Fig. [2] We could also find such critical Richardson 
numbers in the low-Peclet-number regime. 

Moreover, to obtain reliable statistics of the turbulent trans- 
port, the simulation domain must contain several correlation 
length scales of the turbulence and the duration of the simulation 
must be long compared to its correlation time scale. Denoting I 
as the horizontal integral length scale defined by 



' = 271 



(7) 



where E{k\y) is the horizontal spectrum of the turbulent kinetic 
energy averaged over the vertical direction, we find that the ra- 
tio ^/Lh never exceeds 0.16 in our simulations. This indicates 
that there are several large structures in the simulated domain. 
For illustration. Fig. [3] displays 2D snapshots of the temperature 
fluctuations and of the vertical velocity. 

To properly resolve all scales of turbulence, we use a res- 
olution of 128^ X 257 and an aspect ratio Lh/Ly = 7r/4. We 
have verified that in our simulation ^max^ > 1 , where ^^ax is the 
largest wavevector present in the flow, ri - {y^ js)^^^ the Kol- 
mogorov dissipation scale, and s the dissipation rate. This crite- 
rion, which states that all scales down to the dissipation scale are 
resolved, has been found to be adapted for the study of turbulent 
transport in isotropic turbulence simulations ( Gotoh et al. 2002 ). 
In addition, simulations performed at a higher resolution have 
confirmed that our resolution is sufficient. 
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Fig. 3. Snapshots of temperature fluctuations scaled by Ar (left) and 
vertical velocity scaled by At/ (right) in the yOz plane at Pet - 0.34. 
Vertical velocity presents smaller scales than temperature fluctuations; 
both are anti-correlated. 



4, Turbulent transport 

The passive scalar is introduced once a statistical steady state 
is reached. We then consider and compare two complemen- 
tary approaches to determine its vertical turbulent transport: one 
Lagrangian, by following fluid particles; the other Eulerian, by 
solving an advection/difl'usion equation for a concentration field. 

Owing to the stationarity and the spatial homogeneity of the 
turbulence, a fluid particle encounters statistically homegeneous 
conditions as it moves with the flow (as long as it keeps away 
from the upper and lower boundaries). The turbulence is thus 
homogeneous from a Lagrangian point of view. This is an im- 
portant property because it enables us to apply Taylor's turbulent 
transport theory ( [Taylor|192T] ). Accordingly, the transport is dif- 
fusive, and the vertical turbulent difl'usion coefficient reads as 



A 



Jo 



{W(t)w(t + T))dT : 



(8) 



where Tl is the Lagrangian correlation time and () the ensemble 
average over the particles. If z(t) denotes the vertical position of 
a particle, the vertical dispersion is then given by 

(k(0-z(0)]') = 2A^. (9) 

We computed the turbulent diff'usivity using either 
Eq. ([8]) or a linear regression of the mean square displacement 
(k(0 -^(0)]^)- In both cases the main source of error is the 
temporal fluctuations of the averaged turbulent quantities. Then, 
depending on the starting point of the time average, this creates 
a significant dispersion (up to 20%) in the values of A- We fa- 
vor the turbulent diff'usivity obtained with Eq. ([8]) because in 
this case the dispersion is generally lower than using the linear 
regression of the mean square displacement. 

In the Eulerian approach, the concentration field c is gov- 
erned by the equation 

dc 

— + u • Vc = D^Ac, (10) 
ot 

where the molecular diff'usivity is such that Pe^ - 
LbJJ ID^ - 10^. To determine A, we impose a mean concen- 
tration gradient dC/dz and compute 



A 



(11) 



where refers to the concentration deviation from the mean pro- 
file C(z). 
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Table 1. Results of the different runs 



Pee 


Re^ 


Ricv 






P = Ddiui) 


52 


260 


0.124 




6.45 


0.104 


0.90 


150 


1.10 




0.990 


4.01 X 10-2 


0.72 


240 


1.07 




0.773 


9.29 X 10-2 


0.34 


340 


1.27 




0.432 


0.138 


« 1 


335 




0.426 




0.131 



We found that the relative difference between computed 
from the Lagrangian approach (Eq. [8]) and the Eulerian one does 
not exceed 15%. This difference is partly due to molecular diffu- 
sion, which is present in the Eulerian approach but not in the La- 
grangian one. Nevertheless, this effect is limited because /)m/A 
remains always less than 5%. Again, the difference is mostly due 
to temporal fluctuations of the turbulent quantities that generate 
dispersion in the time averages. A way to reduce this disper- 
sion is to average over a longer time. While it is not possible 
in the Lagrangian approach since the displacements of particles 
are limited by the upper and lower boundaries of the domain, a 
much longer integration time can be used in the Eulerian one. 
We have thus used the Eulerian determination of the turbulent 
diffusion coefficient, keeping the error linked to the time average 
to a value lower than 5%. 



5. Results 

We present the results of five simulations corresponding to 
five different values of the turbulent Peclet number, Pe^ = 
52, 0.90, 0.72, 0.34, and « 1 (the last using the SPNA). Table [T] 
displays the turbulent diffusion coefficients, together with other 
relevant physical parameters. 

The most striking feature is that the full computation at 
Pe^ = 0.34 and the SPNA computation give similar results. 
Indeed, at Pe^ = 0.34 the critical Richardson number multi- 
plied by the turbulent Peclet number is very close to the critical 
"Richardson-Peclet" number of the SPNA computation. This 
is strong evidence that at small Peclet number, the steady, sta- 
bly stratified, sheared turbulence is characterized by a critical 
"Richardson-Peclet" number (RiPe{)cY independent of the Peclet 
number. This regime has not been reached in the simulations 
from Pe^ = 52 to Pe^ = 0.7 where RicrPe^ > {RiPee\r- 

Table [T] further shows the parameter p = Di/(u£), which in 
the regime where RiPe^ = (RiPe{)cY, verifies 



Dt = PKRr\RiPef)c 



(12) 



This parameter clearly has the same value for Pe = 0.34 and 
Pe <^ I. Thus, in the small-Peclet-number regime, the turbulent 
diffusion coefficient can be written in the form of (V2\ with a 
constant JS. As a consequence, A is proportional to kRi'^. 

It is obvious that the final expression found for the vertical 
turbulent diffusion coefficient A (see Eq. [12]) has the same form 
as Zahn's prescription given in Eq. ([T]), thus validating in the 
Pe < 0.34 range the hypotheses made to derive it. The present 
numerical simulations also provide a quantitative determination 
of the critical turbulent "Richardson-Peclet" number (RiPe{)cr = 
0.426 and of the coefficient J3 = DJ{u£) ^0.131. This has to 
be compared with the order of magnitude estimates proposed by 
Zahnl ( [T992l ), namely 7?/c < 1/4 andyS = 1/3 ^ 0.333. Never- 
theless, it is the proportionality constant between A and KRi~^, 



i.e. the product of p and {RiPei)^^, which is relevant in practice 
in stellar evolution codes. We find a value of 5.58 x 10-^ for this 
quantity, whereas Zahn's prescription leads to 8.33 x 10-^, which 
is a relative difference of 33%. 

6. Conclusion and discussion 

We thus conclude that the results of our DNS agree with the form 
of Zahn's prescription and provide a first quantitative determina- 
tion of Di beyond Zahn's order of magnitude estimate. In gen- 
eral, the main drawback of numerical simulations of turbulent 
flows is the unrealistically low value of the Reynolds number. 
While it may be interesting to test the dependence of our results 
on the Reynolds number, we note that estimates of effective tur- 
bulent Reynolds numbers associated with the radial diffusion of 
chemicals in stars lead to surprisingly low values Re^ ~ 100 
(Michaud & Zahn 1998), which are comparable to the values 
obtained in the present work. 

Other relevant issues can be addressed numerically. In par- 
ticular, the concentration field is considered here as a passive 
scalar so that the mean gradient has no effect on the flow, 
whereas in reality such a stratification would enhance the ef- 
fect of temperature stratification. This effect is considered dif- 
ferently by |Maeder & Meynet| ( p^96l ) and |Talon & Zahn| ( pW7] ), 
and simulations might help in deciding between them. Another 
important aspect of the models of rotationally-induced mixing is 
the magnitude of the horizontal turbulent transport driven by the 
horizontal differential rotation. It controls the efficiency of the 
meridional circulation transport and the departure from shellular 
rotation (Zahn 1992 ). Again, the turbulent transport can be in- 
vestigated with DNS similar to the present ones by varying the 
angle between the velocity and temperature gradients. 

To study the turbulent transport in less idealized flow config- 
urations would require a better understanding of the shear flows 
induced by the large-scale motions driven by rotation. 
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